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ABSTRACT 

In this paper we propose a new method to locate bed 
boundaries by carrying out a 1-D nonlinear inversion 
of electromagnetic (EM) logging data. We first solve 
for 1-D resistivity structure in which the earth is mod- 
eled using layers of constant thickness. This thickness is 
determined based upon the tool resolution and the de- 
sired resolution from the user. We use general measures 
of data misfit and model structure in the inversion to 
construct piecewise-constant models through the itera- 
tively reweighted least-squares (IRLS) procedure, and 
we minimize a generic global model objective function 
subject to data constraints. The general measure includes 
the traditional l p norm as a special case. The recovered 
piecewise-constant models simulate traditional models 
comprised of a few uniform layers and hence permit easy 
determination of bed boundaries. The generic model ob- 
jective function allows the incorporation of prior geolog- 
ical information and provides measures of the closeness 
to the reference model and the amount of structure in 
the recovered model. We develop an efficient scheme to 
incorporate the IRLS procedure into the search for the 
appropriate regularization parameter at each iteration. 
The iteration? required by the IRLS algorithm arc com- 

r.rc.i.--;'t;: - .:. j'tic ■ip|>jfct!biJi?y of liie «j!*oriibrn in Geijvin- 



INTRODUCTION 

Electromagnetic (EM) logging data are functions of resistiv- 
ity structure and geometry factors such as the length of invasion 
zones and layer boundaries. Inversion techniques can provide 
information about subsurface geological structures. In a typi- 
cal 2-D inversion of resistivity logging data, the earth is mod- 
eled, utilizing the cylindrical coordinate system, as the bore- 



hole, the invaded zone, and the virgin zone. Theoretically, the 
inversion may recover both the geometric parameters (bore- 
hole diameter, invasion length, bed boundaries, etc) and the 
resistivity value for each zone. However, the geometric and 
resistivity parameters have different physical units and affect 
the data in fundamentally different ways. Moreover, the equiv- 
alence between geometric and resistivity parameters results 
in nonunique models, which are difficult to recover simulta- 
neously. One way to reduce the nonuniqueness and cut the 
computational cost of 2-D inversions is to determine the bed 
boundary positions before performing the inversion and then 
to keep the boundaries fixed during the inversion. This method 
has been widely used in the logging industry with success. Most 
existing boundary-picking algorithms use inflection points of 
the input data curves. These algorithms are computationally 
efficient but hi many cases may not meet the required accuracy 
for2-D inversions. In this paper, we propose an alternative way 
to detect bed boundaries through 1-D vertical inversion of re- 
sistivity logging data. In this vertical 1-D inversion the earth is 
modeled using horizontal layers parallel to the surface of the 
earth, and the borehole effect is ignored. 

One-dimensional vertical inversion using an I 2 norm has 
been addressed by several authors (Lin et al., 1984; Dyos, 
1987; Freedman and Minerbo, 1991; Zhang et al., 1994), Re- 
covered models from traditional least-squares solutions are us- 
t:a]J v sw.rtJi representations of the true models. This smeefh- 
:;ci'S, however, dues nui fcc£i;&'t£ nccr.ro le estimation of bed 
h-.y Kit' r: \:: ^:Cir::o:.\ rnj- r.KZV.'Yt of i."-:stfi !T- rohuft 
in the presence oi outliers in the data. Therefore, we investigate 
the solutions of the nonlinear inverse problem using measures 
other than the typical l 2 norm. By judiciously choosing mea- 
sures of data misfit and model structures, we build algorithms 
that are immune to outliers in the data and generate piecewise- 
constant models. 

The use of robust measures of the data misfit in geophysical 
problems has led to significant improvements in EM data pro- 
cessing (Egbert and Booker, 1986; Chave et aL, 1987; Chave 
and Thomson, 1989; Sutarno and VozofX 1991). A few authors 
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have also investigated the application of a non-/ 2 measure of 
model structure in the solution of the nonlinear inverse prob- 
lem. Dosso il990) minimizes the f, norm of the vertical deriva- 
tive of the model to construct piecewise-constant models for 
the 1-D magnetotelluric (MT) problem. He also uses an h mea* 
sure of misfit and solves the linearized inverse problem by using 
linear programming (LP). Oldenburg and Ellis (1993) and Ellis 
et al. (1993^) use an ^ measure to construct blocky models for 
the 2-D MX problem. These authors use LP to obtain their so- 
lutions, Farquharson and Oldenburg (1998) use an iteratively 
reweighted least-squares technique to solve the nonlinear in- 
verse problem of transient EM data, where general measures 
of the model structure and misfit are used. 

In this paper, we develop a new algorithm to determme bed 
boundaries and formation resistivity by solving a 1-D nonlinear 
inverse problem using general measures of model structure and 
data misfit. The traditional l p norm can be considered a subset 
of this general measure. We use horizontal layers with constant 
thickness to model the earth. This thickness is determined by 
the desired resolution from the user and the resolution power 
of the logging tools. The resolution of the recovered bedbound- 
aries is only aa accurate as the chosen layer thickness. The bed 
boundaries aie fixed during the inversion, and the constant but 
unknown resistivity value within each layer is recovered from 
the inversion. A global model objective function is minimized 
subject to data constraints, and the iteratively reweighted least- 
squares (IRLS) technique is used to solve the inverse problem. 
The algorithm is tested on both synthetic and field data. 

METHODOLOGY 

Induction tools consist of individual transmitter- and 
receiver-coil pairs. Responses from each transmitter-receiver 
pair can bd summarized to generate the overall response from 
the tool. In induction logging, a signal generator produces an 
ac current, which is injected into the transmitter coil. The cur- 
rent in the transmitter establishes a primary magnetic field, 
which generates eddy currents in the earth. Tnese eddy cur- 
rents in turn induce a secondary magnetic field in the earth, 
which is a function of physical properties of the earth, and the 
receiver coils measure the voltages induced by the secondary 
magnetic field. The tool response is further converted to ap- 
parent resistivity defined as (Moran and Kunz, 1962) 



(■a = 



;\*T H t : 



where NT Is the number of turns of the ith transmitter; N ; 
is the number of turns of the ;th receiver; V iS is the induced 
voltage received by the ith receiver from the yth transmit- 
ter with smgle turns; and N T and N R are the numbers of the 
transmitters and receivers, respectively. The value L is given by 
V T?\ Nt N il L ih where L '> fa ^ ^ separation of the 
ith transimtter and the j th receiver. Trie final data are given m 
ohm-met en . 

Since our intention is to detect the bed boundary, we are 
less concerned about the accuracy of the resistivity. Therefore, 
in the forwird modeling, we neglect the borehole and inva- 
sion effects. This reduces the forward modeling to a 1-D prob- 
lem which is solved following the work of Anderson (1986) 
andKenrtedy et al. (1986) for an arbitrarily oriented magnetic 



dipole in a horizontal-layered earth. To locate the bed bound- 
aries accurately, we model the earth with a relatively large num- 
ber of layers of constant thickness. This thickness is determined 
according to the tool resolution or the desired resolution from 
the user. 

A GENERAL MEASURE FOR VECTORS 

A typical geophysical inverse problem involves minimizing 
a model objective function and/or a data misfit objective func- 
tion. These model and data objective functions are obtained 
by taking the norms of the model structure vector and the data 
misfit vector, respectively. Different norms can be used as a 
measure of a vector, whether that vector defines the misfit be- 
tween the observed and predicted data or describes a model. 

Consider a vector x with elements XjJ = 1. 2 lv\ The 

traditional l p norm of x is given by 

where p > 1. When p = 2 this expression is the h norm of the 
vector. In this paper we propose a new measure 

where p > 0 and * \xj I is some positive number. We explain 
how to choose € later using numerical examples. Zhang et al. 
(1998) have proven, from a statistical point of view, that p in 
equation (3) must be > zero. When p > h the above general 
measure reduces to the one proposed by Ekblom (1973, 1987) 
and used by Farquharson and Oldenburg (1997). As * becomes 
small and p > 1, the measure in equation (3) tends to the l p 
norm in equation (2). . 

The parameter p in the general measure controls weights on 
the elements of a vector. Ib illustrate this we create a vector 
x = (xx x 2 x N ), whose elements range from 1<T* to 18. We 
then compute {x)+€ z Y* for all components of this vector 
with p « 2, 1, 0.5, and 0 (Figure 1) . Here we set e to a negligible 
value of 1Q" 6 . As p decreases, the weights on small vaiuesof 
x f increase and the weights for large values of xj decrease. The 




Fig. 1. Comparison of general measures with different values 
of p. Trie value of € is 1(T*. 
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general measure is the sum of each (xl + € 2 Y n \ so the smaller 
p is, tJie more important the small components in ||x||£. Setting 
p == 0 leads to a uniform weighting, which equals 1 for all values 
of Xj, Hie resultant general measure equals the dimension of 
x. In nhis extreme case, all elements of x contribute the same 
toward the general measure ||x||£. 

Until now, no attempt has been made in the literature to use 
equation (2) as a measure of misfit or model structure with 
0 S P < 1, where the l p norm is not defined. From a.numerical 
point Df view, when p < 1 the derivative of ||xP has a singularity 
at xj == 0, 1 < ; < N . Our general measure can be considered a 
perturbation and extension of the traditional i p norm, arid the 
singularity at xj 0. 1 < j ' < N is removed with the introduced 
small positive constant e. Below, we use this general measure 
to develop our inversion algorithm. 

INVERSION ALGORITHM 

In the nonlinear inverse problem we are solving, the model 
m(r) is the spatial distribution of resistivity of the earth, dis- 
cretizcd as 

" . - v« • 

«n(r) = 2^m y ^(r), s (4) 

where N is the number of the model parameters; r is the po^ 
sition vector; and $j is the basis functions. In this paper, we . 
chcosii y/f j as boxcar functions. Suppose the forward modejiag; 
problem, which calculates the response from a given models 
expressed as . 

where m= (mi , m 2 ffijv) r ; d pfe is the data vector generated •: 

from this model; and f is the nonlinear mapping that connects 
the model space to the data space. The inverse problem bpijh 
sisto of finding, the model from the collection that reproduces 
the observed data to a desired misfit level, which is Ae most 
plausible given any prior geological and geophysical ^forma- 
tion. ./ 

Since this is a nonlinear inverse problem, we use lineariza- 
tion arid iterate for the solution. The objective function to be 
minimized is 



(6) 



<t> = fa + P~ 4>d- 
The model objective function is given by 

*■= a*. +(1 -«)«,. (7) 

where fi, = ||W,(ni - ni«,)|]?. and = ||W x (m - m^Jllf . The 
data objective function is given by 

*4 | [W d (d° ta - d"" 1 ) - J(m"- 1 )5m] [|f - 4>?. (8) 

where d"" 1 is the predicted data from m* -1 , the model frpm 
the previous iteration, and J is the sensitivity matrix, 



J u (m- 1 ) = 



dntj 



(9) 



where in^ is the reference model that can be used to incor- 
porate a priori geological information into the solution. The 
parameter is the target misfit level at each iteration, which 
is reduced gradually — typically by a factor between 2 and l(h- 
to prevent the violation of the linearization. The value 0 is the 



regglarization parameter. In equation (7), the two terms on 
the right side provide (a) the measure of the closeness of the 
recovered model to the reference model and (b) the amount of 
structure in. the reconstructed model, respectively. The value 
.W, is the weighting matrix associated with the minimum-norm 

: solution. In this paper, we assume the reference model is of 
equal importance over the entire spatial extent and thus choose 
the identity matrix as W x . Since, a minimum-norm solution 
may possess unwanted, detailed structures when noise exists 
in the date; a;jEb$t-order, finite-difference operator W x is used 

. to rru'riimi^ jtiie; amount of structure presented in the recov- 
er^ }nib^b^e instant a controls the balance between the 
minimum-norm and minimum-structure solutions. When a = 0 

4 and p tstZi^is. is a least-squares solution that has the minimum 
ntpdcl horrrt Trie smoothing operator W x can be introduced 
by^ertj^O ^ <1, with containing statistical information 

;^^U^'^^.iif^ia.:3Che noise in the data is assumed to be uncor- 
Telate^,^ W^js diagonal If the noise is Gaussian with zero 
: . meander! tiie diagonal elements of W rf are the reciprocals of 

; the,standtircl deviations, s 

•^^^ffiP^^^cQ^^^ Action in equation (7), we set 
•; |^||^^|^^ijer entiate;. <X> with respect to the elements of 

llA^*ija^^^^t6 v te resulting equations to zero. The system of 
. • J;jpqu<iil^nsAw b.q sdjyedis given by 





where x lt u it and v, are the ith elements of the corres- 
ponding vectors \V,(m - m nt ). W^re - m,*), and Wrffd 1 * 5 - 
d(m!'" 1 ) - J(ra"~ )Sm] t respectively. The detailed derivation of 
equation (10) is given in the Appendix. We explain how to 
choose values of € later using numerical examples. When p 
tends to zero, the system in equation (10) reduces to that with- 
out regularization because and R, are both zero matrices, 
'^^upe :^ e : diagonal matrices R„ R„ and so that 
.-K*>^i^gft^V.«rf Rjfirf^Rrf, then equation (10) 




(14) 



^MI^^^Ai R„ and R,, are all diagonal, 
■■^^iS^^^^^^^P 1110 ^ easily. 



at 



. . : , 
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Equation (14) is very similar to the system of equations that 
would be obtained if all measures in the objective function 
were sum-of-squarcs measures. Hie iterations required by the 
IRLS procedure have been absorbed into the iteration associ- 
ated with the linearization of the nonlinear problem. At each 
iteration associated with the linearization of the problem, the 
system of equations given in equation (14) is constructed and 
solved usingthe re weighting matrices R,, R x , and R rf computed 
from the model obtained from the previous iteration. 

The regularizauon parameter /3 is determined in the follow- 
ing manner. Be fore the final desired misfit is achieved, we solve 
equation (14) at 2 to 4 logarithmically spaced values for fi in a 
relatively small interval, usually between 0.001 and 100, at each 
iteration. We do not update matrices R, , R x , and R^ in the line 
search. Once the final target misfit for the inversion is reached, 
the regularizatlon parameter fi is fixed and further iterations 
are carried out to allow the JRLS procedure to converge. The 
iteration is termbiated when the change in the model norm 
reaches a predetermined value, usually set to 1%. 

NUMERICAL EXAMPLES 

In the following;, we demonstrate the benefits of using generic 
measures of model structure and data misfit by applying the 
inversion to synthetic dual-phase induction log (DPIL) data 
generated over a two-layer resistivity model. The DPIL tool is 
a conventional induction logging tool that operates at 20 kHz. 
It consists of several transmitter-receiver pairs and provides 
medium and deep resistivity curves. Since the medium curve 
has higher vertical resolution, we use it for boundary picking. 
The true model and the synthetic medium curve for the model 
are shown in Hgure 2. The data were generated for a single 
frequency of 20 IlHz at a station interval of 0.1 m. For conve- 
nience, we refer 1 o the simulated data from the true models in 
synthetic examples as the observed data and the data computed 
from the models recovered from the inversion as the predicted 
data. 

In the inversion, the earth was parameterized into 100 layers 
of 0.1-ro constant thickness. This thickness is determined from 
experiments and is about the limit of the resolution from the 
inversion of DPIL data. Such a parameterization not only al- 
lows us to determine the bed boundary locations accurately but 




also provides us with a fine vertical variation of the resistivity 
structure. A starting model of a 5-ohm-ro half-space was used 
for the following two examples. Such an initial model is easy 
to build. However, better starting models (if available) can be 
used to reduce the number of iterations in the inversion. The 
parameter a was chosen to be 0.002 based on our experience. 

EM data are usually contaminated with noise. Two common 
types of noise are Gaussian noise and outliers. Our experience 
indicates that, while a small amount of Gaussian noise has little 
effect on inversions using general measures, outliers can affect 
inversions significantly. To illustrate the benefits of using ro- 
bust measures of misfit, we perturbed two data samples in the 
synthetic data in Figure 2. The data values at 3 m and 8.2 ra 
were changed from 10.56 to 20 ohm-m and 5.25 to 10 ohm-m, 
respectively. 

We first inverted the perturbed data using l 2 measures as 
the model structure and data misfit The best-fit solution was 
reached after eight iterations, which corresponds to a 6% rms 
error. The recovered model from the inversion (Figure 3a) 
has been affected severely by the two outliers in the data, and 
the overall structure is very rough. The shallower resistive zone 
has shifted upward, and an artificial resistive layer has been 
created between 8 and 8.5 m. The predicted data in Figure 3b 
are trying to fit these two outliers. This example highlights 




4 5 
Depth (m) 

Fig. 2. The true model (the solid line) and the synthetic data 
to be used in the inversion (the crosses). 



4 6 
Depth (m) 

Fig 3 Results from the inversion of the perturbed data using 
the> 2 norm, (a) The model. The solid line denotes the recovered 
model and the dashed line denotes the true model (b) The 
predicted (solid line) and observed (plus signs) data. 
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the vulnerability of an l 2 measure of misfit to outliers in the 
data. 

We then test our algorithm using p = 1 on the same data set 
Fust, we determine the values of e for computing <p 3 , 4> Xi and 
fa in equations (7) and (8). These values of 6 should be small 
enough so the final values of fa, fa, and fa will not differ 
significantly from their values with f =0 yet large enough to 
make the IRLS procedure converge. Before carrying out the 
inversion, we do not know the final values of fa, fa, and fa 
necessary for determining values of €. In a borehole environ- 
ment, however, the observed data, which are given in apparent 
resistivity/conductivity, arc a reasonable estimation of the tine 
resistivity structure. Thus, we can use the data in Figure 2 to 
calculate fa, fa> and fa as functions of e (Figure 4). From these 
three curves, we choose values of € that are small enough so 
that fa , (f> x , and fa differ only in the third significant figure from 
then values with € =» 0 yet are large enough to make the IRLS 
procedure converge. The values of <f for the model structure 
and fcr data misfit are set to 0.01 and 10 -5 , respectively. 

After 20 iterations, the best rms misfit of 65% is reached. 
Figure 5 shows the model and data from the inversion. The re- 
cover^ model is an excellent representation of the true model. 
The two outliers in the observed data are, in effect, ignored by 
the predicted data. This indicates that the fit to the spiky data 
obtained using essentially an h measure of misfit is more ap- 
propriate than that obtained using an l 2 measure. This example 
also indicates that the use of an /j measure of model structure 
can lead to more blocky models from the inversion. 

Tb further test our algorithm, we apply it to a 2-D syn- 
thetic data set generated from a modified Oklahoma bench- 
mark model shown in Figure 6. This model has a large range of 
layer thicknesses and resistivity values with conductive inva- 
sion. The synthetic data are generated at a sampling rate of four 
sample s per foot for a DPEL tool. Data from 80 to 270 ft (24 to 
82 m) are used in the inversion, and the subsurface is parame- 
terized with0.5-ft (0.15-m) thick layers. A 1-ohm-m half-space 
is used as both the starting model and the reference model in 
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FIG. 4. The variation with e of general measure (with /? = 1) 
for & (curve 1), fa (curve 2), and fa (curve 3) computed using 
me data in Figure 2. r & 
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F*0. 5. Inversion of the perturbed data using a robust measure 
for the misfit and model structure, (a) Model, (b) Data. 
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Fro. 6. Oklahoma benchmark model, (a) L xo (length of inva- 
sion), (b) R xo (resistivity of invasion), (c) /?, (resistivity of for- 
mation). 
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» inversion. The parameters a, p. and « far the misfit and 
,odel structure are set to 0.002, 0.5, and If 5 -"^^ 
ruall threshold value of 0.001 is used for the change of model 
orm in the IRLS procedure to obtain maximum blockmess 
„ the recovered model. After 20 iterations, the inversion con- 
•erges to the final model in Figure 7a. The use of a general norm 
1S the measure of model structure in the inversion enables the 
•Lovery of the p.ecwise^onstaut model, which m turn allows 
fasy deletion of the bed boundaries. The estimated reastm- 
S are also close to the true resistivity values. Because of the 
resolution limitation of the deep curve, the inversion does not 
U,We the thin lay,** around 240 ft (73 m). ™'™™"^™ n 
sistivity model is an average of the invaded and true bnnabon 
resistivity and hen«differs from the t™e model. The predicted 
and observed data in Figure 7b agree within 5% rms error. 

] 1ELD DATA EXAMPLE 

Asafield data example, we applied ^ e P ro P^^°iJS 
te the dual-phase induction data acquired from the Abu Dhab. 
region in the Middle East. In the inversion, we used the deep 
data from 8465 to 3600 ft (2580 to 2620 m), acquired at a sam- 
pling rate of four samples per foot. The subsurface was param- 
eterized using 0-5-ft (0.15-m) layers and a starting 
of 1 ohm-m. Inversion was performed using p = 0.5, « = omu. 
and an « value of 10- for both misfit and model structure. 
In the line search, misfit evaluations were conducted at four 
logarithmically equispaced points between 0.001 and 100/ The 
inversion was performed on the deep induction curve. There 
wasanexcellentdatamatchbetween the field andtoe^nAetic 
data obtained from the inversion (Figure 8b) to about 1% rms 
error. The final model obtained from the inversion in Figure 8a 
seemingly resolved most of the major formation layers. 

CONCLUSIONS 

The proper locntion of bed boundaries and good prior infor- 
mation about the resistivities of the formations, among various 

other factors, play a key role in the suceess of a 2-D resisW 
inversion scheme. Conventionally, bed boundaries are located 
based on the slope variations in the borehole EM dau. This 
may result in unrealistic bed boundaries, especially in regions 
of strong shoulder-bed effects. Inversion schemes can be used 
to estimate the tied boundary locations and to obtain accurate 
a priori informat.on about resistivity for 2-D inversions. 
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Fir 7 i^True ( dashed line) and recovered (solid line) models. 
fb) Obsirv^i fi&gri predicted (solid line) data from the 
inversion of a 2-D deep induction curve. 
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r rp^hIm from the inversion of an Abu Dhabi deep- 
Sducuon'urve (e^vered model, (b) Observed (dashed 
line) and predicted (solid hne) data. 



In this paper, we have developed and implemented a model- 
„orm-baL P in^ 

fit and model structure, in which a generic global model ob- 
jective function is minimized. A large number of Am layers 
lere used to accommodate arbitrary variations of the geo- 
logical structure An iterative reweighted ^ t 
durewasusedinmeinversiontogenerateapiecewjs^nst^^ 

model that allowed easy determination of bed boundaries The 
results obtained from the synthetic and field data ^aqta 
proved the effectiveness of the proposed scheme. The mver- 
Ln of synthetic log data showed mat ourmveniona^onti^ 

L recover the bed boundaries and resistivity of *e Jorma- 
tions effectively, even in the presence o l™^*^??™* 
the algorithm on a field data example from the MiddteEast 
demonstrated the ability of the proposed inversion ^ algor dim 
to recover the bed boundaries and resistivities £ J* «™j 
tion from field data Though the proposed algondun * not a 
replacement for 2-D inversion, it provides excellent a .priori 
nformation for the 2-D inversion. Also, in regions where the 
subsurface is one dimensional (no invasion present), the a go- 
S works effectively in estimating the resistivity variation 
of the subsurface. 
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APPENDIX A 
DERIVATION OF EQUATION (10) 

EOS to! ° f ^^terminequation (7). Identical reasoning leads to 

-n-w,(m- ). Then the general measure of x is 

A „ itoi — = *T RxW,(m»-i -Mm - m Kl ), 

where Uie elements of the diagonal matrix R, are given by 

= i (gww., - , j g^.^tr, (A . 5) 



(A-l) 

tend"! o = m " ' + 6m differeDtiatin 8 the above equation 



3||W rf [d Qb * - d"~' - J(m"-i)3 m ] ||f 
dSm 

= J r WjR <f W (< [Ji n ,4-d obs -d"- 1 ]. (A-6) 
9.5m" = ^W** K" 1 + *"» ~ rare,), (A-2) where the elements of the diagonal matrix R„ are given by 

where R, is a diagonal matrix whose elements are given by P |[g<^) y [<» - /yCra"- 1 )*™;]^ 

The summary of equations (A2), (A4), and <A6) and reorga- 
nization yields equation (10). 



f T AT -|2 IP/2-1 

'{[EWvK +A . ( A . 
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